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^ ; ABSTRACT 

^ \ Among the most explored directions in the study of dense stellar systems is the 

\ investigation of the effects of the retention of supernova remnants, especially that 

^\ of the massive stellar remnant black holes (BH), in star clusters. By virtue of their 

O I eventual high central concentration, these stellar mass BHs potentially invoke a 

I wide variety of physical phenomena, the most important ones being emission of 

1 gravitational waves (GW), formation of X-ray binaries and modification of the 
^ ' dynamical evolution of the cluster. Here we propose, for the first time, that rapid 

>• ' removal of stars from the outer parts of a cluster by the strong tidal field in the 

2 ' inner region of our Galaxy can unveil its BH sub-cluster, which appears like a 
^ ■ star cluster that is gravitationally bound by an invisible mass. We study the 
^ ' formation and properties of such systems through direct N-body computations 

I ' and estimate that they can be present in significant numbers in the inner region of 

^ ■ the Milky Way. We call such objects "dark star clusters" (DSCs) as they appear 

• dimmer than normal star clusters of similar mass and they comprise a predicted, 

■ new class of entities. The finding of DSCs will robustly cross-check BH-retention; 

5^ ■ they will not only constrain the uncertain natal kicks of BHs, thereby the widely- 

debated theoretical models of BH-formation, but will also pin-point star clusters 
as potential sites for GW emission for forthcoming ground-based detectors such 
as the "Advanced LIGO". Finally, we also discuss the relevance of DSCs for the 
nature of IRS 13E. 

Subject headings: stellar dynamics — methods: N-body simulations — open 
clusters and associations: individual (IRS 13E) — Galaxy: center — gravitational 
waves — black hole physics 
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Introduction 



Compact remnants of massive stars in star clusters, which are neutron stars (NS) and 
black holes (BH), form a dynamically interesting sub-population due to their tendency of 
segregating towards the cluster's center and augmenting their population density therein. In 
this respect, the BHs are special in that they undergo a "runaway" mass segregation. These 
remnant BHs are typically several 10s of Mq heavy, enough to form a Spitzer-unstable sub- 
system, provided a significant number of them are retained i n their paren t cluster. Due to 
this instability (also called the mass-stratification instability, ISpitzerlll987| ). the continually 
sinking BHs cannot come to an energy equipartition with the local surrounding stars and 
finally end-up in a central, highly concentrated sub-cluster made purely of BHs , which is self- 



gravi t ating and dyiiamically nearly isolated fro m the rest of the stellar cluster (IMerritt et al. 
20041 : iMackev et allboosi iBaneriee et aPboiol ). 



Such a dense environment of BHs is dynamically very active due to the formation of 
BH-BH binaries vi a 3-body enc ounters (iHeggie &: Hutll2003l ) and their hardening by super- 
elastic encounters (jHeggielll975l ) with their surrounding BHs. Studies of the dynamics of pure 
BH sub-clusters using Monte-Carlo and direct N-body integration methods indicate that the 
dynamical BH-BH merger events they generate are likely to contribute a significant gravita- 
tional wave (G W) detection rate to the future "Advanced LIGO" 



observatories ( 


Porteeies Zwart & McMillan 


2000 




BenacQuista 
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Moodv & Sie;urdsson 


2OO9I: 


Banerjee et al. 


2OI0I; 


Downing: et al. 


2011) 



Such studies show 

that a BH sub-cluster is typically self-depleted in a fe w Gyr due to the super-elastic dy- 



nami cal encounters and the resulting escape of the BHs (lO'Leary et al.ll2006l : iBanerjee et al. 



2OIOI ). Th e energy extracted from the tight BH-BH binaries heats- up and expands the clus- 
ter's core ( IMerritt et al.l 120041 : iMackey et al.ll2008[ ). which can be detectable by the future 
optical missions such as the "GAIA" mission. Furthermore, the BHs can be important for 
dynamically formed BH X-ray sources due to their encounters with the surrounding stars 
( llvanova et al.l 12010). X-ray observations have indicated the prese nce of BH X-ray binary 
candidates in GCs (iMaccarone et al.l 120071 : iBrassington et al.ll2010l ). The presence of a BH 
sub-cluster within a star cluster, therefore, has the potential to give rise to a plethora of 
physical phenomena, all of which have significance to upcoming prime missions such as the 
GAIA, AdLIGO and the present and future X-ray missions. 

Is it possible to obtain any direct observational signature of the presence of a BH sub- 
cluster within a star cluster? We predict here, for the first time, that within a few kpc from 
the Galactic center, rapid tidal stripping of star clusters by the strong tidal field can expose 
its BH sub-cluster. This would happen when the timescale of the preferential removal of 
stars from the outer regions of the cluster is shorter or comparable to the encounter-driven 



- 3 - 



self- depletion timescale of its central BH sub-cluster (see above). Such a dissolved phase 
of the cluster would consist of a few stars orbiting around a cluster of BHs and would 
observationally appear as a highly super-virial star cluster with a large mass-to-light-ratio. 

As we discuss here, a number of direct N-body computations of model star clusters 
indeed support the formation of such systems. These objects comprise a predicted, new 
class of compact stellar populations which we name "dark star clusters". The importance 
of dark star clusters (hereafter DSC) is twofold: on one hand, if they are found to exist, 
then they guarantee that star clusters are potential sites for GW emission and formation of 
BH X-r ay binaries and on th e other hand, they naturally constrain the uncertain natal kicks 
of BHs (IWillems et al.ll2005[ ). as DSCs can form only if a significant number of BHs retain 
in the cluster following their proge nitor supernoyae. This, in turn, restricts the theoretical 
models of core collapse supernovae (jjanka et al.ll2007| ). 



2. Computations 



We compute the evolution of model star clusters subjected to the Galactic tidal field 
using the direct N-body integration m ethod. For ou r purposes, we use the state-of-the-art 
N-body integration code "NB0DY6" (lAarsethl 12003), which, apart from utilizing a highly 
sophisticated numerical integration scheme (IMakino fc Aarsthl 119921 ). also follows the evo- 
lution of the individual stars until their r emnant phases, using an analytic but well-tested 
stellar evolution recipe (IHurley et al.ll2000l ) . A unique feature of NB0DY6 is its use of highly 



accurate regul arization methods in resolving close encounters (IKustaanheimo fc Stiefellll965 



Aarsethll2003l ). Furthermore, the code exploits the remarkable hardware-accelerated com- 
puting capacity of Graphical Processing Units (GPUs) in integrating the centers of masses. 
NB0DY6 currently incorporates general-relativistic effects only through analytic prescrip- 
tions of GW energy-loss. 



2.1. Dark star clusters 



We follow the evolution of initial Plummer clusters (jKroupall2008[ ) of single stars, having 
masses between lO^M© < Mcz(O) < 7.5 x lO^'M© and half-mass radii between 1.0 pc < 
^h{0) ^ 3.5 pc. All the clusters initially consist of zero-age-main-sequence (ZAMS) stars 
with their masses m chosen from the canonical initial mass function (IMF; lKroupall200ll ) 
p{m) oc m°, where a = — 1.3 for O.OTMq < m < O.5M0 and a = —2.3 (Salpeter index) for 
m > O.5M0. Their metallicities are chosen to be solar, as suitable for our Galaxy's disk. We 
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assume for now that all the supernova remnants [i.e., BHs and NSs) receive low natal kicks 
in general such that they remain bound to their parent clusters at their formation. Such 
models follow circular orbit around a point mass of = 2 x 10^*^ Mq representing the Milky 
Way bulge. 

Fig. [1] (top panel) shows an example of the evolution of the virial coefficient for one of 
our computed model clusters with initially iV(0) = 6.5 x 10^ stars and r/i(0) = 3.5 pc, located 
at Rg = 2.0 kpc Galactocentric distance. The orange curve shows the time-evolution of the 
virial coefficient, Q, for all the members bound to the cluster including the BHs and the NSs, 
which mostly remains constant at Q ~ 0.5, as it should be for the quasi-sta tic relaxation of 
a self-gravitating system through two-body encounters (iHeggie fc Hutll2003l ). 



The green curve in Fig. [T] (top) shows the time-evolution of the virial coefficient taking 
into account only the luminous objects, i.e., the nuclear burning stars and the white dwarfs 
(hereafter WD), which are those an observer sees. The BHs formed have masses ~ IOMq and 
the NSs are typically of ~ 2Mq, which are significantly more massive than the majority of the 
remaining luminous members and hence are segregated to the cluster's center. As the lower 
mass luminous stars in the outer regions are stripped by the external field, the gravitational 
potential of these central invisible remnants becomes increasingly important. Therefore, 
the kinetic energy of the luminous sub-system increasingly exceeds the corresponding self- 
equilibrium (or quasi-static) value as the constituents perceive a potential that becomes 
increasingly deeper than their self-potential. This correspondingly raises their exclusive virial 
coefficient, Q^, above 0.5 as in Fig. [T] (top panel, green curve). The cluster thus evolves to 
a DSC state (see Sec. I2.1.ip . Note that while reaches a very high value, the system as a 
whole remains bound since Q < 1 throughout (except at the final dissolved state, not shown 
in the figure). 

Notably, NSs form a few factors more often than the BHs and contribute significantly in 
elevating Q^. However, the NSs being lighter than the BHs, their sub-population occupies a 
more extended zone in the cluster's center. This causes them to get stripped earlier than the 
BHs {c.f. Fig.lll bottom panel) so that in the late evolutionary phase, typically when Q^, > 1, 
it is mostly the BH population that contributes to the augmented Q^,, thereby determining 
the hfetime of the DSC phase (see Sec. I2.1.ip . Nevertheless, at large enough distances 
from the Galactic center, where slower tidal stripping cau ses the DSC phase to appear 
later than the self-depletion of the BH sub-cluster (see Sec. [H Sigurdsso n fc Hernoui st*^1993l 
Merritt et all I2OO4I lO' Learv et all bood : iMoodv fc Sigurdsson! 12009^ .Downing et al. .20101 : 
Banerjee et al.ll2010l ). the NSs constitute the primary dark component of the corresponding 
DSC state (see Sec. I2.1.ip . This self-depletion process, although operative for both the 
NS and the BH sub-populations, is more efficient for the latter as it is significantly more 
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concentrated. 



2.1.1. Galactic population of DSCs 

Fig. [2] (top) shows the expected increasing trend in the hfetime of the DSCs with initial 
cluster mass M^iifS) {Rq = 2 kpc, rh{0) = 3.5 pc). The DSC phase can be defined when the 
cluster appears unbound, i.e., > 1.0 or when it appears significantly super- virial, which 
we take when Q* > 0.75, and denote the corresponding lifetimes by tdscq* with = 1.0 
and 0.75 respectively. Fig. [2] (bottom) shows tdsc,i-o against Rq {Mci{0) = 3 x IO^Mq, 
r/i(0) = 3.5 pc). Tdscq, ^^^^ increases with increasing Rq since the DSC takes a longer time 
to get depleted in a weaker external field. Notably, for Rq > 4 kpc, the DSC state becomes 
NS-dominated. Beyond Rg ^ 5.5 kpc, it takes more than the age of our Galaxy 10 Gyr) 
for the representative cluster to evolve to its DSC phase {i.e., Q^, > 0.75). 

It then follows that for the present-day Galaxy the DSCs are formed from clusters with 
Mci{0) > 10*^ M0 within Rg ^ 5 kpc, to be taken as representative numbers, and have 
lifetimes tdsci.o ~ 150 Myr and tdsco.ts ~ 250 Myr as conservative estimates. Although 
these estimates are based on a point-mass tidal field, note that the DSCs' progenitor clusters 
form in the Galactic disk and orbit on nearly circular paths (in the equatorial plane) and 
hence would experience the same external field with an axisymmetric disk-like distribution 
of the same mass. The tdsc.q* — Rg dependence can however be moderately modified as the 
clusters see a mass increasing moderately with Rg- 

To estimate the Galactic population of DSCs, we take the average star cluster formation 



rate (hereafter CFR) over the entire Galactic disk [Rg < 10 kpc) to be O .IGM^^ yr ^ (ILarsen 



2OO9I ) following a Schechter initi al cluster mas s function ( lSchechterlll976l ) over the mass range 
IOOM0 < Mci{0) < 2 X 10^ Mq (jLarsenll2009l ) which is assumed to remain invariant over the 



last few hundred Myr. This implies ~ 0.5 clusters form per Myr which have properties 
(Mc;(0) > IO^Mq and Rg < 5 kpc) that must have them evolve to DSCs if a sufficiently 
high fraction of BHs are retained. Assuming a steady state conversion to the DSC phase 
with lifetime tdsci.o ~ 150 Myr (see above), the expected number of > 1 DSCs in the 
Galaxy within Rg < 5 kpc is A^dsc.i.o ~ 75. For DSCs with Q^: > 0.75, the corresponding 
number is A^dsc,o.75 ~ 125. Hence, a significant number of DSCs can be expected in the 
inner Galactic zone. 
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2.2. Can IRS 13E be a DSC? 



There has been recent concern with the widely debated IRS 13E (hereafter IRS13E); 
an extremely compact stellar association of a few young, massive stars at a close projection 
to the Galactic ce nter that apparently survives the extreme tidal field by being bound by 
an invisible mass ( iMaillard et al.l 120041 ) . While a ~ 1300Mq intermediate mass black hole 



(IMBH) was widely believed to be this invisible mass ( jPortegies Zwart fc McMillan! 12002 



Maillard et al.l |2004J : iPortegies Zwart et al.ll2006i). thi s possibility has recently been ruled 



out with a significant confidence by iFritz et al.l (120101 ) through their newer proper motion 
measurements, leaving the nature of the invisible component of IRS13E currently ambiguous. 
This status of IRS13E prompts us to consider whether its dark component can be an ensemble 
of stellar-mass BHs instead. To that end, we perform preliminary N-body calculations to 
determine whether a DSC configuration resembling IRS13E can be a possible fate of a star 
cluster very close to the Galactic center. 

These computed clusters follow circular orbits around the Galactic SMBH (a central 
mass of ~ 5.4 x lO^M© which is a combination of the mass of the SMBH and that of the 
nuclear star cluster within Rq < 1 pc) within Rg < few pc, where s t ar formation has been 
shown to possibly lead to a top-heavy stellar IMF (|Bartko et al.ll2010l : iNayakshin fc Sunyaev 



20051 : iNavakshin et al.l 12006 



Bonnell fc Ricd |2008|) . Since we are primarily interested in 



the final state, we begin the N-body calculations from an evolved phase of the cluster, for 
computational ease. Therefore, we initiate the computations with Plummer clusters made 
of stars which are pre-evolved until tq ^ 3.5 Myr age. This age is slightly earlier than when 
the most massive star (of ~ I5OM0) evolves to a BH. All the stars are thus still on their 
main sequence. 



Observations indicate a very flat IMF 



0.45 ± 0.3; see iBartko et al.l 120101 ) for 



stars close to the Galactic center and also a dearth of low-mass stars. The latter is evident 



from a declining density of B-stars away f rom the central SMBH (|Bartko et al.ll2010l ) and a 
significant lack of coronal X-ray emission ( iNayakshin fc Sunyaevll2005l ) from the SgrA* field. 
It is currently unclear from the literature what would be the lower-mass limit of such an 
IMF, which, in turn, determines the a of the mass function (MF) and its lower limit Mi at 
cluster-age tq. In this preliminary study, we simply take the total number of stars in the 
cluster A^, a and Mi at age tq as free parameters. The upper MF limit is chosen to be the 
canonical I5OM0 and the stellar metallicity is solar. These limits are the ZAMS values and 
they are appropriately reduced during the tq ~ 3.5 Myr pre-evolution by the NB0DY6's 
built-in stellar evolution prescription (see above). Like the computations in Sec. 12. H we 
retain all the BHs (of ~ lOM© as obtained from within NB0DY6) in the cluster. At the 
beginning of the computations, the models are taken compact to a similar extent as IRS13E 
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( iMaillard et al.ll2004j ). with half-mass radii between 0.02 pc < rh{0) < 0.04 pc and they orbit 



within 2 pc < Rq < 4 pc. 

We find that to reach a IRS13E-hke state, i.e., a state where typically ~ 10 luminous, 
young stars are very tightly bound to a cluster of ~ 100 BHs, it takes ^ 6.5 x 10^ and a 
rather extreme lower cut off of M; ^ 35 Mf^ fo r MF index a ^ (at age tq), this a being close 
to the upper limit of the iBartko et al.l ( 120101 ) index. A steeper MF requires even higher Mi. 



As a demonstration, it can be seen in Fig. |3] (top panel) that this cluster eventually evolves to 
a configuration consisting of ~ 130 BHs, comprising a ~ I3OOM0 dark component, orbited 
by < 10 young stars. Although the half-mass radius of the system expands at the beginning 
of the evolution, primarily driven by the mass loss through the massive stars' winds and 
their supernovae, it finally collapses to ~ 0.02 pc, as shown in Fig. |3] (middle panel), by 
the time the system arrives at the above configuration. The luminous members in this state 
include 1-2 0-stars (those still on the main sequence) and a few helium and Wolf-Rayet 
stars, thereby being of similar variety as observed in IRS13E. The final state of the system, 



there fore, resembles IRS13E in terms of compactness and stellar content ( iMaillard et al. 



20041 ). Fig. [3] (lower panel) shows that the luminous sub-system becomes highly super- virial. 



i.e., a DSC, as the IRS13E-like state occurs. 

The calculations presented in this section are preliminary which suggest an interme- 
diate state of a cluster, close to the Galactic center, that would evolve to an IRS13E-like 
configuration. It remains an open question whether any reasonable initial cluster evolves 
to such an intermediate phase, which depends on the low-mass limit of the Galactic-central 
IMF, its index and the initial mass and compactness of the cluster, at a given Rq. A scan 
over these parameters, beginning with much larger and more compact clusters, is necessary 
to determine such possible initial configuration(s), which is much more compute-intensive. 
Although the above Mi appears too high for a lower cut-off, its progenitor cluster would have 
a smaller limit taking into account the rapid tidal stripping over the 3.5 Myr pre-evolution. 
Tidal stripping also implies a steeper initial index, in better agreement with the observed 
one. 

Although we over-simplify by effectively ignoring the dynamical evolution of the cluster 
during the stellar pre-evolution, the dynamical history is not instrumental in determining the 
occurrence of a IRS13E-type state which happens merely due to the competition between 
stellar evolution and tidal dissolution of the cluster. It is enough that the cluster remains 
bound by the time most of the stars become BHs quenching the wind mass loss; the system 
would then core-collapse to a compact configuration irrespective of its history. This is in 
contrast with the "classical" DSCs discussed in Sec. 12.11 whose formation depends crucially 
on mass segregation. 
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A potential drawback of this "in-situ" model is that it is likely to significantly over- 
produce the number of massive stars seen in the Galactic center. 

Given these drawbacks, the calculations in this section are only to suggest that IRS13E 
might be a DSC but are by no means conclusive. 



Discussion 



Our calculations (Sec. 12. ip signify that a gravitationally bound star cluster naturally 
evolves to an apparent super-virial state, while remaining bound as a whole, as a consequence 
of the interplay between the dynamics and the evolution of stars, provided a significant 
number of stellar remnants survive in the bound system after their formation via supernovae. 
The existence of such dark star clusters is a first-time prediction and serves as an excellent 
cross-check of the retention of supernova remnants in star clusters, the effects of which 
is widely explored. The very presence of DSCs would kill two birds with o ne stone by 



havin g consequences on the widely debated theoretical stellar collapse models (jjanka et al. 



20071 ) due to the implied direct constraints on natal kicks, and by securing s tar clusters as 
potential sources for the forthcomi i ig "Advanced LIGO" ( iHarry et al.ll2010l ) GW detector 
( jPortegies Zwart fc McMillan! l2000t iBanerjee et al.l 120101 ). Given that the DSCs must have 
a significant population in the inner region of our Galaxy, it is tempting to conduct a survey 
of intermediate-aged (see below) stellar assemblies to identify them. 

Dark star clusters can be observationally distinguished from actually dissolving clusters 
(which are also super virial and hence expanding and are generally young) through their 
compact sizes in spite of their intermediate ages, as shown in Fig. HJ That the apparent 
super-virial state of a DSC is not due to an IMBH, would be indicated by the absence of a 



centr al cusp in its velocity-dispersion profile, which would be there otherwise (iNoyola et al. 
2008h . 



The present studies justify DSCs as predicted, new type of compact stellar populations. 
An immediate improvement over this study is to explore any effects of varying orbital eccen- 
tricity and different initial profile types {e.g., using King instead of Plummer profiles). Also, 
the effect of a varying star (and hence cluster) formation rate over the Galactic disc needs 
to be incorporated to determine the predicted Galactic DSC population, as an improvement 
over our assumption of a uniform average cluster formation. Another development would be 
to consider a disk-like mass distribution instead of a central point-mass (but see Sec. I2.1.ip . 



The preliminary computations in Sec. 12.21 suggest that the dark component of IRS13E 
can be an ensemble of stellar mass BHs and IRS13E may therefore perhaps be a DSC. 
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However, this conclusion should, for now, be taken as being suggestive rather than conclusive 
due to the drawbacks discussed (Sec. 12.21) . We postpone a more detailed and self-consistent 
study on this issue in a future paper. 
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Fig. 1. — Top: Evolution of the virial coefficient of a Plummer cluster with A^(0) = 6.5 x 10^, 
r/i(0) = 3.5 pc and Rq = 2.0 kpc. The lower (orange) curve represents the virial coefficient, 
Q, for the whole cluster, i.e., including the luminous stars, NSs and BHs bound to the 
cluster. The upper (green) curve represents the virial coefficient, Q*, of the luminous objects 
exclusively. The monotonic growth of implies the approach towards an apparent highly 
super-virial state — the "dark star cluster" phase. Bottom: The corresponding evolution 
of the numbers of the luminous stars, A^*, the NSs, A^ns and of the BHs, A^bh- 
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Fig. 2. — Top: The lifetimes, tdscq*; DSCs with increasing initial Plummer cluster 

mass Mci(O) {Rq = 2 kpc, rh{0) = 3.5 pc). Bottom: tdscq. (Q* > 1) as a function of Rq 
for a Plummer cluster with Mci{0) = 3 x lO^M© and r/i(0) = 3.5 pc. 
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Fig. 3. — The evolution of the numbers of the luminous members, N^, and of the BHs, A^bh, 
for a model star cluster computation, leading to an IRS13E solution. The direct N-body 
computation starts with a Plummer cluster of pre-evolved stars having A^(0) = 6.5 x 10^, 
'"/i(O) = 0.025 pc and Rg = 4 pc, where the stars are taken from a a = MF within 
35M0 < M < I5OM0 ZAMS mass (see Sec. \2.2\\ . The evolution of the Lagrangian radii 
of the system is shown in the middle panel. The final state of the cluster resembles that 
of IRS13E. The evolution of the viral coefficient Q, for the whole system, and Q^,, for the 
luminous sub-system is shown in the bottom panel. 
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Fig. 4. — The evolution of the Lagrange radii (top) and the core radius (bottom) of the 
example cluster of Fig. [1] The half-mass radius (black curve, top panel) remains nearly 
constant throughout and the core also remains compact following an initial brief contraction 
phase. 



